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Stability Analysis in Magnetic Resonance 
Elastography II 


Heiko GirnperleirU Alden Waters* 


Abstract 

We consider the inverse problem of finding unknown elastic param¬ 
eters from internal measurements of displacement fields for tissues. In 
the sequel to [4], we use pseudodifferential methods for the problem of 
recovering the shear modulus for Stokes systems from internal data. We 
prove stability estimates in d = 2, 3 with reduced regularity on the es¬ 
timates and show that the presence of a finite dimensional kernel can 
be removed. This implies the convergence of the Landweber numerical 
iteration scheme. We also show that these hypotheses are natural for 
experimental use in constructing shear modulus distributions. 
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1 Introduction 

This article uses pseudodifferential methods to sharpen recent stability estimates 
for an inverse problem of Magnetic Resonance Elastography (MRE), with short 
proofs. The new proofs allow to analyse practical numerical reconstruction 
methods. 

In Magnetic Resonance Elastography, internal measurements of time-harmonic 
displacement fields offer the possibility of a highly resolved reconstruction of 
shear modulus distributions. It is motivated by the detection of cancerous 
anomalies in their early stages m See [5J [6J 1 EH HI EH M for relevant 
related works. 
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The reduced regularity assumptions of the stability estimates in this article 
prove relevant for numerical reconstruction schemes. We show how they directly 
relate to classical results for overdetermined elliptic boundary problems and 
their analysis using pseudodifferential operators. The analysis is based on a 
Stokes system as in g] . 

Let f1 denote a simply-connected bounded domain in where d = 2, 3 with 
C°°-boundary dVl. We consider the following boundary value problem for the 
elasticity equations 


T V(A(a;)V • u\) + uj 2 u\(x) + 2V • p(x)V s u\(x) = 0 in fl, 
(_ u\(x) = F(x) on dfl , 

where 


( 1 . 1 ) 


2 V s u\ = Vu\ + (VuxY 

and Vua is the matrix (djU\,i)f j = \, with u\^ the i-th component of u\. The 
Lame coefficients (respectively, the shear and the compressional modulus) A, /i £ 
C 1 (fi) satisfy 

A > A min = min{A(a;) : x £ fi} > 0, (1-2) 

p > /iram = min{^(a;) : x £ Cl} > 0. (1-3) 

If F € ff 1//2 (dfi) and (11.21) and (11.31) are satisfied, there exists a unique solu¬ 
tion u\ £ Lf 1 (fl) d to (11.11) even for A, fi £ In particular, V s ma £ L 2 (Vl) d . 

Higher regularity V s u\ £ H 4 (Q) d holds under the additional assumption that 
A,/x£C 4 (H), FeH 9 / 2 {dn) d . 

In g] it was shown that for 2/i max < 3A m i n the solution u\(x) is approxi¬ 
mated in IL 1 (H) ci norm up to 0( A -1 / 2 ) by the solution to the Stokes problem 

a i 2 u(x) + 2V • /r(x)V s u{x) + Vp(:r) =0 in O, 

V • u(x) = 0 in H, 

< (1.4) 

u( x) = F{x) on dfl, 

f p(x ) dx = 0. 

If we examine the solutions in 2 and 3 dimensions, we find we can recon¬ 
struct a displacement of the shear modulus fi in a stable way from the difference 
of the solutions. Hybrid modalities involve exciting the system with more than 
one wave or modality. In 3 d it was already shown in g] using elliptic regular¬ 
ity theory that hybrid modalities are necessary for reconstruction of the shear 
modulus. This necessity is in contrast to dimension 2. The goal of this article 
is to show how pseudodifferential operators and the theory of elliptic boundary 
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problems leads to simple proofs of stronger stability estimates for the inverse 
problem. The estimates reduce the necessary order of regularity and allow the 
analysis of practical numerical reconstruction methods. 

In Section [2] we introduce basic notions of elliptic boundary problems and 
the Shapiro-Lopatinskii ellipticity condition. A basic Fredholm theorem is re¬ 
called in Theorem [T] In sections [3] and [I] we apply Theorem [I] to conclude basic 
stability estimates (Theorems [3] and [2]) in the presence of a finite dimensional 
kernel. Refined stability estimates in L , based on stronger assumptions, are 
the content of Section 5. We also show that locally it is possible to remove 
the presence of the finite dimensional kernel under certain hypotheses on the 
system, Corollaries [3] and [I] These results imply convergence of the numerical 
Landweber iteration scheme, which is discussed in Section [8] We also show 
that the assumptions which are made on the symbols are natural for numerical 
experiments in Section [7] 


Notation: In this paper we use the Einstein summation convention. For two 
matrices A and B , the inner product is denoted by 


A : B — aijbji, 

and we write \A\ 2 = A : A. For vector-valued functions 


f{x) = f 2 {x),f d (x)) : ft ->• , 


the Hilbert space H m (iT ) d , m G N is defined as the completion of the space 
C^°(£l) d with respect to the norm 


E/ 

N=i n 


(|V */(*)| 2 + l/ 0 )l 2 ) dx , 


where we write V 1 = d 11 .. . d ld for i = (*i, ..., id) for the higher-order deriva¬ 
tive. 

The outward unit normal to the boundary df2 is denoted by n. If /i £ C 1 (f2), 
we define the conormal derivative 

2 d -£ = M (a0 (V/ + (V/)‘) n . 

We say that a linear operator T : C'“(K. d ) —t C°°(K d ) is of order < a if 

\\Tu\\ m <C\\u\\ m+a , «€HR d ), 

for every met. We say ord(T) = a if T is of order < a, but not of order < a' 
for a' < a. 
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2 Reduction of the Regularity for Stability Es¬ 
timates 

This section recalls some aspects of the elliptic theory for boundary problems in 
low-order Sobolev spaces. These will be applied to obtain regularity estimates 
for the inverse problem which are suitable for numerical calculations. 

We outline a direct approach to Fredholm estimates for a class of boundary 
problems relevant to the reconstruction problem. Our presentation will follow 
the treatment in Eskin [7], Ch. 52, to account for the low-order Sobolev spaces. 

After straightening out the boundary, it is sufficient to consider the boundary 
problem in the half-space 

R” = {(x',x n ) : x n > 0, x' £ M”' 1 }. 

We consider A(x, D) as an r x r system of differential operators of order to, 
with principal part Aq(x, D). We say than an operator A(x, D) is elliptic if the 
determinant of the principal symbol detAo(x,£) defines the symbol of a scalar 
elliptic operator of order mr. In other words, we have detA 0 (x,£) ^ 0 for all 
£ ^ 0. We consider the initial boundary value problem 

A(x, D)u = f in 0 

Bj(x',D)u\dn = gj 1 < j < m+ = to/2 

We let B{x ', D) be the m + x r matrix with rows Bj. The system of boundary 
conditions can then be written in the following form 

B(x',D)u\an = g, 

with g = (< 71 ,.., g m+ )- We can associate a family A 0 (x' , £' x ,, 0, D n ) and B 0 (x' , £' x ,, 0, D n ) 
which are differential operators on depending on (x' ,£' x ,) £ T 0 * (<9f7). 

The Shapiro-Lopatinskii condition can be stated as the unique solvability 
of the differential equation in L 2 ( R+). In order to verify this condition, we 
translate it into a non-vanishing condition on a determinant associated to the 
symbols of the differential operators. 

If Aq(x, 0, + 1 ) = 1 , we obtain a factorization Aq(x, C, fn) = A+{X, £n) A- (X, £', Cn), 

with 

m_|_ 

A±(x,C^n) = ~ cr j ± ( a: ^'))- 

i=l 

For any N > s an integer, we construct operators with symbols Rp (x, ^ n ) 

and degi?+ = —m+ — p for |^'| > 1 and operators with symbols A~(x, 
with degA“ = m+ — p with |/'| > 1 and 0 < p < N such that 

( N \ N 

J2 R${x, D~,D n )\=Yl ( X ’ D ' D ") + T m+-N- 1 . 

P—0 / p=0 
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The operators T mp _|_ at_i are of lower order: 


ll^mpH—iV-l M llo >s < C IMIm + -l,s-iV 
Because A is a differential operator, the product 


Vs. 


A(x, D)Rp (x, D', D n ) 


is of the following form 


m 


3=0 


where 


Bo p (x, £) = A 0 (x,€)R+(x,£',£ n ), 

Bjp(x,0 = J2 1 9 & D k x R+ + Aj(x,Z)R+ , j> 1, 


r —0 \k\=j—r 


k\ d£ k 


A{x,£) = Aq(x,£) + ^2Aj(x,£) ordAj = m-j. 

3=1 


If we take the product of A and J2p=o Bp an d collect symbols of the same degree 
of homogeneity in (£',£„) with |£'| > 1 we obtain 

A 0 (x, £)R£ (x, £„) = A~ (x, £„), 

Ao{x,£)R+(x,£',€n) + T p (x,€',£ n ) = A~(x,£',£ n ), p > 1, 
with T p depending on Rf for 0 < j < p — 1, 


T P = 


j+r=p 


B, 


and degTp = m + — p for |£'| > 1. 

We now define the operator IT as 

n’D(x',0,t',t n ) = ^ j Dd£ n , 

7 + 

with 7 + being a contour which encloses the poles of D in the upper half plane 
and let 


bjk(x\0 = Il'B jk (x',0,£',£ n ). 

Here Bj & is the symbol of the composition of Bj and R + ^ k _ 1 . 
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The Shapiro-Lopatinskii condition in its algebraic form then says that 

de% fco (^£% m + = i^O V^'UVO. (2.1) 

Here bjko(x',f') denotes the prinicipal part of bjk- 

We let A(x, D) be elliptic in Q and the Shapiro-Lopatinskii condition (12.11) 
be satisfied. We let 


s> max (toj + 1/2), 

1 <j<m + 

with rrij = deg Bj. Under these conditions, we then have the following theorem: 

Theorem 1. The boundary value problem defines a Fredholm operator from 
H s (fl) to 

m_|_ 

H s —m X 1/2 - 

1=1 

There exists a constant C such that 

IML.n — C ^11/ll-^.n + 5Z • 

A similar theorem can be found in the earlier work [8], but we use the 
formulation of [7] for the ease of computation of the symbol classes. 

3 Stability in Dimension 2 

We prove the following analogue of Theorem 3 in :|]. 

Theorem 2. Let fl C R 2 a smooth and bounded domain, (u\,p\), ( U 2 ,P 2 ) be 
two solutions to equation m3 with coefficients pi, ^2 E C 4 (U), respectively. 
Assume that p\ = p 2 on dfl and 

|V s ui(a:)| ^ 0, xen. (3.1) 

Let e G (0,1). Then there exists a finite dimensional subspace K C L 2 (Ll) and 
a constant C > 0, depending on H/T 2 IIc 5 / 2 + e (f 2 )> such that 

ll/ll — p 211!/2+e,f2 — ^ IK — U2 lli/2+e,n > (3-2) 

provided p± — p 2 -L K. 

In |5 under the same hypotheses we proved 

ll/ii -M 2 || 4 ,n <C\\ui -u 2 \ t 6>n 

The difference between this theorem and Theorem 4 in [4] is the norms on the 
estimates in 15.11 and the shorter proof. The reduction in regularity is more 
suitable for numerical experiments as discussed in Section [8] 
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Proof. For the proof, we may assume without loss of generality that d Xl u\(x ) ^ 
0. We apply the operator (d Xl , —d X2 ) to (11.41) in order to eliminate the pressure 
term. For g = gi — /X 2 , we derive 

(d xi , —d X2 ) ■ (2V • gV s ui) = g. (3.3) 


By definition, g satisfies 


IMIi.O — ^ ll u l M 2||;_|_2,n 

for all l < 0. The principal symbol of the linear operator in (13.31) is given by 

A)(a:,0 = 2|C| 2 d Xl u\(x). 

In order to verify the Shapiro-Lopatinskii condition ED, we compute the zeros 
of the polynomial in terms of £i. Provided d Xl u(x) ^ 0, the root in the upper 
half plane can be computed as 


Ci = ±*£2- 


It follows that for j = 0,1 


bj 0 - 


2ni J z — 1 
7+ 


t dz ^ 0. 


We then apply Theorem [Q 


□ 


4 Stability in Dimension 3 

We prove the following theorem, which is the analogue of Theorem 3 in . 

Theorem 3. Let (ui,pi) and (ui,pi) be solutions to ED with different bound¬ 
ary conditions. In other words we set F(x) = Fi(x) and F(x) = Fi(x) in ED 
for the respective solutions but they share g = g\. We assume that there ex¬ 
ists a positive constant C independent of (x,f) £ T*Ll, where T*Q denotes the 
cotangent space, such that 

c|£| 4 > |(VV0z)£) X £| 2 + KV-firtoO X £| 2 > i|£| 4 . (4.1) 

Let (u 2 ,P 2 ) and ( v, 2 ,p 2 ) be solutions to the Stokes system ED with g = g 2 
and F{x) = F\(x) and F{x) = F\{x), respectively. Let e £ (0,1). Assume 
that gi,g 2 £ C 7 / 2+e (fl) and g\ = g 2 on dfl. Then there exist a constant C, 
depending on ||M2||c 7 / 2 + e (n) an d a finite dimensional subspace K of H l ^ 2+t (Ll) 
such that 

WTI ~ 4^211 i/2+e,0 — C (jl M l — U 2||3/2+e,f2 ~HI^1 ~ ^21 l3/2+e,o) , (4-2) 

provided that (gi — gf) _L K. 
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The main difference from Theorem 3 in [5] is the reduction of regularity 
on the norms. However there is a slightly stronger assumption on the symbol 
classes, in [4] the upper bound was not necessary. 

We begin as previously by eliminating the pressure terms from the Stokes 
systems. We consider the equations for i = 1, 2 which are 

V • fj,i\7 s Ui + u 2 Ui + Vpi = 0 

When we take the cross product of both sides this yields the equation 
V x V • fiiV s Ui + w 2 V xuj = 0 

If we set p, = pi — H 2 , w = u\ — U 2 and subtract the first equation from the 
second equation, we obtain 

A Ul H = V x V • hS7 s u\ = g (4.3) 

with 


g = —V x [V • /x 2 V s wj] — w 2 V x w. 

Clearly we then have that there is a constant C which depends on ||At 2 |lc 2 + I (? 2 ) 
norm such that 


IIfl'lli,o — C IIHIj+3,n 

for all l. We view the identity (|4.3I) as an over-determined second-order partial 
differential equation with principal symbol 

(V- Ul (z)0 x £ 

which is not elliptic. We augment the operator with a second set of measure¬ 
ments: 


A U2 /j, = V x V • hV s U 2 = g. 

We set 

A * = (Ad Ai) B * = (Trace on 90,Trace on dil), 

and we analyse the new system 

AM = (AM,AM) = (g, g, 0,0). (4.4) 

The condition (14.111 ensures that the system is now injectively elliptic, whence 
AA is elliptic and satisfies the assumptions of Section [2] If (AA)~ is a 
parametrix of A A, a left parametrix of A is given by (AA)~A- 

Therefore, if the augmented system defines an elliptic boundary value prob¬ 
lem, we conclude: 

IMkn < C'dl^lU—2,^2 + ||s|| s - 2 ,n) < C'dlui - U2||«+i,n + ||«i — U2||s+i,n) ■ 

To show the Shapiro-Lopatinskii condition (12. II) . we argue as follows. Ac¬ 
cording to Hypothesis (14.11) . the polynomial |(V s ui(x)£) x£| 2 + |(V s ui(x)£) x£| 2 
is elliptic and therefore, if r = Vsl+W ^ 0, two of its roots each lie in the 
upper resp. lower half-plane. As the roots are homogeneous in r, it is easy to 
show that det (bjk o(x', £ , ))™ A )' ^ 0. We then apply Theorem|T] 



5 Reduction of Regularity in Weighted Sobolev 
Classes 

To further reduce the order in the stability estimates of Theorems [2] and [3l we 
recall the Fredholm properties of a Shapiro-Lopatinskii elliptic boundary value 
problem, A = (A, B) with ord(T) = m. See QT| for more detailed results. 

Fix a boundary defining function p £ C°°(Q). For s £ No we let 

0 s = {u : p^D a u £ L 2 Va, |a| < s}. 

We set 0 _s = (0)*. Elements in 0” s are of the form 

/= T, D a (p^f a ) for some f a £ L 2 . 

\oc\<s 


Let 


D° A {n) = {u€L 2 :Au£ 0~ s (fl)} 

endowed with the graph norm of A. The basic Fredholm theorem then reads: 
Theorem 4. A defines a Fredholm operator from D A (fil) to 

m .|_ — 1 

0- m (fi) x H H s - m i- 1/2 (T). 
j =0 

This allows us to conclude in dimension 2: 

Corollary 1. Under the hypotheses of Theorem 0 there exists a finite dimen¬ 
sional subspace K C L 2 (fl) and a constant C > 0 depending on pi and tt such 
that 


Wui - T2\\ L 2 {Q ) < C\\p 2 {u 1 - M 2 )|| L2(n) (5.1) 


provided pi — P 2 A K. 

In dimension 3 we have: 

Corollary 2. Under the hypotheses of Theorem 0 and there exists a constant 
C, depending on and O and a finite dimensional subspace K of L 2 (il) 

such that 

IImi -HI L2(n) < c(||p^ 2 («i -w 2 )|| ffl(n) + ||p _2 (wi - w 2 )|| H i (n) ), 

provided that (pi — P 2 ) -L K. 
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6 Local Injectivity Results 

Recall that we say the principal symbol is non-characterstic at the origin if 
p(x ,£) is such that p(0,£) ^ 0 with £ = (0,0,0,1). We set r = We 
assume the following three hypotheses, given some e > 0 

1. The symbol p(x,£',t) = 0 has at most simple real zeros and at most 
double complex zeros in terms of r. 

2. If ti,T 2 are distinct zeros of p(x,£',r) = 0 then this implies |ti — T 2 I > e. 

3. For any non-real zero r of p(x,£',r), we have |3 t| > e. 

We recall the following uniqueness result, based on a refined Carleman esti¬ 
mate [ISl- 

Theorem 5. Assume that the plane x n = 0 is non-characteristic at the origin 
and the Hypotheses (1-3) above hold. Let l > 0 and w € H l (R d ) a solution of 
Pw = 0 in a neighborhood of the origin which vanishes identically for x n < 0. 
Then w = 0 in a full neighborhood of the origin. 

We prove: 

Corollary 3. We assume in d = 2 that u(x) 7 ^ 0. If p = 0 in some 
neighborhood of the hyper-plane x n = 0, and u\ — U 2 = 0, then p = 0 in a 
neighborhood of the origin. 

Proof. We recall that the principal symbol is 

m 2 diy. 

Setting t = £ 1 we are looking for roots of the equation r 2 = — £f. These 
solutions r = ±*£2 clearly satisfy the necessary hypotheses. We must also check 
that the polynomial associated to Aq(x,£) is non-characteristic. The hypothesis 
is satisfied under the assumption df. 1 u(0,x') 7 ^ 0 , which we already assumed in 
order that the operator be elliptic. □ 

We also prove: 

Corollary 4. We assume in d = 3 that the condition holds. If p = 0 in 
some neighborhood of the hyper-plane x n = 0, and u\ — 112 = 0 and u± — U 2 = 0, 
then p = 0 in a neighborhood of the origin. 

Proof. We recall the principal symbol 

IA,MIH(v s «iOxOI + l(v 5 %Ox£|. 

For Aq(x, () to be noncharacteristic, we compute the symbol. We need only that 
one of the off-diagonal terms in the matrices defining the symmetric gradients 
V s ui and V s Mi be nonzero at the origin. The roots of the polynomial satisfying 
conditions (i-iii) are determined by the condition EH), as discussed in Section 

m □ 
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7 Examples Illustrating the Boundary Condi¬ 
tions 

We would like to show that there are natural means to enforce the somewhat 
unusual looking hypotheses on the symbols given by m in d = 3. 

Lemma 1. It is possible to chose initial data F and F to the system (II.dl) such 
that the condition (HUD holds a.e. 

Proof. We let u\ be a solution to the system 113). Here we note that V s «i 
satisfies a closely related elliptic problem 

w 2 V s iti — V S (V • /iiV s Mi) = V s Vp (7.1) 

with corresponding boundary conditions. The boundary data V s Mi|an is deter¬ 
mined by F through a Dirichlet-to-Neumann map. We chose F to be a random 
field on dfl such that V s Mi|ao is a random field of Wigner matrices. In other 
words, the entries of V s rti|an are distributed according to a Gaussian law and 
trV s wi|afi = 0. Because u\ is a solution to E3, it is a given that trV s iti = 0 
in ft. Because the solution operator R which is associated to m is a lin¬ 
ear pseudo-differential boundary operator, X/ S u(x) = R(X7 s u |an) still obeys a 
Gaussian law for every x £ 12. It follows as a result of m that 

P(v : 3A £ R + , Av = Xv, Bv = Xv) 

for A , B with A ^ B Wigner matrices is small. In other words, £ is an eigen¬ 
vector of V s u and V s u with small probability provided V s ui|an and V s uigQ 
are chosen as Wigner matrices. This necessarily implies the condition 14.11 with 
a high probability. The natural condition in dimension d = 2 is that V s u > 0, 
which can be enforced by similar means. □ 

8 Stability of the Landweber Iteration Scheme 

We consider as in @J the reconstruction of the true shear modulus distribution 
Ptr- The treatment follows the book [5] closely. We want to reconstruct the /i 
from u m which is the measured displacement field. We introduce the functional 



n 
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where we have that u is the solution to El which minimizes J[p\ when p. is 
varied. Fixing p we consider a solution v to 


2V ■ pV s v + ui 2 v + X7p = (u — u m ) in fl, 

V ■ v = 0 in 12, 

< 

v = 0 on dfl, 

J pdx = 0. 

We then compute the Frechet derivative DJ[p] of J as 

< DJ{p\, Sp >= J 5pX7 s v : S7 s udx. 
n 

This identifies the functional DJ[p\ withV s u : V s u. 

With the gradient descent method, the numerical minimization of J amounts 
to the following. After an initial guess po, we update it with the following 
scheme: 


p n+1 (x) = p n (x) - CrDJ[p n \(x) X G fl, 71 > 0, (8.1) 

with a being the step size. This procedure is outlined in [2]. 

Following [J], the mapping T 


T : [i u, 


is such that 


dj\iA = (DF\fi]y(r[n] - J'KD, 

where the superscript * denotes the adjoint. 

The resulting optimal control scheme (18.111 is controlled by a Landweber 
iteration scheme given by 

Hn+l{x) = Unix) - Cr(DT[p\)* (^[p] - T[fltr])(x) X G fl, 71 > 0. 

From 1 Appendix A], the following convergence result in L 2 (f2) for the 
Landweber (or equivalently the optimal control) scheme holds. 

Theorem 6. Let d = 2. Assume that the assumptions of Theorem\^are satisfied 
and K is trivial. If, for sufficiently small e 0 and some e G (0,1), 

HMD — l x tr||^fi/2+e(Q) < £0i 

then 


|\hn ~ A*tr | |^i/ 2 +€(q) — > 0 as 71 —>■ OO. 
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The modification between this theorem and the result in [2] is the reduction 
of the necessary regularity to practically relevant values. The discrepancy func¬ 
tional in d = 3 must be modified to include both sets of data as in [2] . We can 
then conclude that: 

Theorem 7. Let d = 3. Assume that the assumptions of Theorem\^are satisfied 
and that K is trivial. If, for sufficiently small e 0 and some e £ (0,1), 

11 Mo — M*llfji/2+e(n) < e 0i 

then we have 


IIMn — Mtr||#i/ 2 +e(Q) —> 0 as oo. 
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